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Abstract This work presents a technique for particle size generation and 
placement in arbitrary closed domains. Its main application is the simulation 
of granular media described by disks. Particle size generation is based on the 
statistical analysis of granulomctric curves which are used as empirical cumu- 
lative distribution functions to sample from mixtures of uniform distributions. 
The desired porosity is attained by selecting a certain number of particles, 
and their placement is performed by a stochastic point process. We present 
an application analyzing different types of sand and clay, where we model the 
grain size with the gamma, lognormal, WeibuU and hyperbolic distributions. 
The parameters from the resulting best fit are used to generate samples from 
the theoretical distribution, which are used for filling a finite-size area with 
non-overlapping disks deployed by a Simple Sequential Inhibition stochastic 
point process. Such filled areas are relevant as plausible inputs for assessing 
Discrete Element Method and similar techniques. 

Keywords granular media • simulation • particulate systems • particle 
generation • point processess 

1 Introduction 

The Discrete Element Method (DEM) has become a powerful tool in the 
numerical simulation of engineering problems involving discontinuous media. 
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Among the problems where the method has been apphed, we may cite fragmen- 
tation, fracture, impact, and coUision phenomena, in addition to those directly 
related to soil modeling as part of studies of geomcchanical problems [19,20, 
47,48]. 

DEMs are based on media discretization into finite sets of particles. More 
often than not, this discretization must satisfy some requirements. In order 
to achieve significant results, a good representation of the studied media is 
needed and, therefore, the particle generation must meet the granulometric 
curves of soil and should be as close as possible to the specified porosity. A 
general procedure would be then comprised of two stages, namely, (i) particle 
size generation, and (ii) particle placement. 

The measurement of particle size is an activity that is common to diverse 
disciplines: archaeology, fuel technology, medicine, geology, sedimentology, civil 
engineering, pharmacology, etc. The common objective is to determine the 
overall size distribution of a collection of particles [24]. Particle-size distri- 
bution (PSD) is fundamental for characterizing construction materials, soil 
mechanics, soil physics, sediment-flux in rivers, and others. In soil science, 
it is typically presented as the percentage of the total mass of soil occupied 
by a given size fraction. Determination of the soil PSD is not a trivial task 
because of the heterogeneity of the shape and density of particles [22]. Clas- 
sical techniques used to determine the PSD, like sieving and sedimentation, 
are point-wise and, thus, require an interpolation to obtain the complete PSD 
curve. The transformation of discrete points into continuous functions can be 
made by mathematical models [21], where the normal, lognormal and WeibuU 
distributions are especially prevalent [24]. 

The PSD in soils is frequently assumed to be approximately lognormal [13, 
43], but there are soils which present bimodal PSDs [45]. A bimodal lognormal 
Gaussian distribution is presented in [42] for the characterization of various soil 
samples. It consists of a weighted sum of two distributions: primary minerals 
(sand and silt), and secondary minerals (clay), each described by a Gaussian 
law. 

Besides the lognormal distributions, other models have been proposed, in- 
cluding one based on the water retention curve [26], the Gompertz model [37], 
the fragmentation model [12] and estimating the PSD from limited soil texture 
data [44]. Furthermore, the modelling of PSD by means of the fractal mass 
distribution is presented in [35]. Combining some well-founded theoretical re- 
sults from fractal geometry, the model allows to simulate the PSD of a given 
soil and its characterization by means of the entropy dimension. 

Although statisticians have ocasionally examined particular problems with 
the estimation of PSD [25], it was not until the work of [5] that a coherent 
statistical approach was formulated. That paper introduced the log-hyperbolic 
distribution as a suitable model for particle sizes. However, the considerable 
computational difficult involved in fitting this model has prevented it from 
achieving widespread use [24]. 

The comparison of mathematical models for fitting PSD curves in soil sci- 
ence has been performed in a few works [14,21,27,41]. On the studies of [27], 
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seven parametric models were tested - the majority lognormal models and 
growing curves type - including a unimodal and six bimodal models. They 
compare five lognormal models with one, two and three parameters [14], the 
Gompertz model with four parameters [37] and the Fredlund model with four 
parameters [26] using 1387 Corean soil samples. Four comparison techniques 
were considered to define the best model: the coefficient of determination (i?^), 
the F statistic, the Cp statistic of Mallows [34] and the Akaike's information 
criterion (AIC). They conclude that the Fredlund model presented the best 
fit for the majority of the soils studied, with increasing performance with the 
increase of clay content. Furthermore, they showed that texture could affect 
the performance of the PSD models. This work represents an important con- 
tribution to the model comparison for fitting granulometric data, but another 
models, potencially adaptable to this finality were not studied. Moreover, the 
authors compared models with different number of parameters; this could be 
a problem because it can benefit those more degrees of freedom [21]. 

The work done by [21] tested and compared fourteen different models with 
feasibility to fit the cumulative PSD curve based on four measured points. 
The parameter used to compare the models was the sum of the square errors 
between the measured and calculated values. They concluded that the most 
recommendable models to fit PSD curves were: Skaggs model [44], WeibuU 
model [4C] and Morgan model [36], all of them with three parameters. 

In this work we transform data from granulometric curves into unevenly 
spaced histograms, and then sample from this empirical distribution. The sam- 
pled data is explained by means of four distributions, namely the hyperbolic, 
gamma, log-normal and Weibull laws, and the best fit is associated to the 
original granulometric curve. We then obtain an arbitrary number of radii 
sampling from the chosen distribution, and we place these particles in an arbi- 
trary closed region of the plane or the space, using a stochastic point process. 
The number of particles to be used is iteratively determined in order to satisfy 
a predefined porosity. 

Many procedures have been presented for particle placement with well- 
known tractable shapes, especially disks and balls (spheres) for the two or 
three- dimensional cases, respectively. These procedures are generally divided 
into two distinct groups: geometric and dynamic algorithms. 

Dynamic algorithms use mechanical procedures that reproduce external 
actions to achieve an initial set of particles. A typical approach of these al- 
gorithms consists of creating a regular array of particles and them randomly 
disorder it through the action of gravitational forces [23]. Other approaches 
require the use of DEM simulation procedures in order to achieve the initial 
configuration. 

These dynamic approaches have many advantages, including the possibil- 
ity of filling any arbitrary-geometrical domain with elements of predefined, 
irregular sizes and, still, controlling the desired porosity. An overview of typi- 
cal dynamic techniques can be found in [18,23]. However, in most cases, these 
dynamic algorithms require checking of contacts between particles to achieve 



4 



A.C. Frery et al. 



the desired result, making them significantly slower in terms of computational 
cost. 

Another category of techniques is based on geometric algorithms. These 
algorithms are characterized by the use of purely geometric concepts for gen- 
erating particles in a given domain, without the need of dynamic simulations 
to represent their movement and interactions. This often reduces significantly 
the processing time. The literature of granular media offers several techniques 
for random generation of disks or spheres with different diameters. One of the 
main strategies consists on establishing the position and the size of particles 
from random numbers. If there is overlap, a new random location is performed 
with the same fixed size [28,32]. 

Other geometric techniques are based on particle generation from an initial 
triangular or tetrahedrons elements mesh, and placing disks or spheres in the 
interior of these elements [18]. Several other geometric strategies can be found 
in the literature, including those that use the concept of boundary contraction, 
where the particle locations are calculated from the previous inclusion of other 
particles with a pre-established diameter [2,31,33]. 

The class of Random Sequential Adsorption - RSA processes describes the 
deposition of particles. Variations include the dimension and other properties 
of the substrate, the shape of the particles and how they interact among them 
and with the media [38]. An important characteristic of our work is that we 
use random particle sizes. Quoting Cadilhe et al [15]: 

Versions of the continuum RSA model with deposition of mixtures of 
particles of different sizes and shapes have only been studied for the 
simplest cases of two particle sizes and identical shapes. [. . . ] In spite 
of all the research work thus far, the field remains widely open to new 
research efforts. 

Dynamic algorithms grant the simulation of samples which are in geostatic 
balance with spatial homogeneity, and either specified granulometry or prosify, 
at the expense of high computational cost. Geometric algorithms are much 
faster than the former, but they do not grant geostatic balance. The latter 
are, thus, more adequate in, for instance, Monte Carlo techniques. Moreover, 
traditional geometric algorithms have difficulties in reaching answers that meet 
the specified granulometric curves. 

This work presents a geometric technique for particle size generation and 
placement in arbitrary closed domains. Its main application is the simulation 
of granular media described by disks which can be extended to spheres. Parti- 
cle size generation is based on the statistical analysis of granulometric curves. 
These curves are treated as empirical cumulative distribution functions, his- 
tograms are derived from them, and each bin is used as the density of a uniform 
random variable. Samples from these uniform laws are then obtained. The de- 
sired porosity is attained by selecting a certain number of particles, and their 
placement is performed by a stochastic point process. This proposal grants 
both porosity and granulometric properties with spatial homogeneity, i.e, the 
system is in global equilibrium as in fully saturated porous media (see Fig- 
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ure 1). On the one hand, by construction, the results adhere to the statistical 
properties of the observed data. On the other hand, the porosity we measured 
is an underestimate of the true property; nevertheless, the method accepts and 
produces any value of porosity, so it is up to th 
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Fig. 1 Fully saturated porous media (From [16]) 

The proposed strategy is comprised of the following steps: 

1. Particle size generation: 

(a) A granulometric curve and the porosity are used as input (section 2.1). 

(b) A histogram is derived from this curve, and samples are drawn from 
the empirical distribution provided by the histogram (section 2.2). 

(c) The samples are used to estimate the parameters of a set of distributions 
(section 2.3). 

(d) The estimated distribution that best fits the original data is adopted 
(section 2.4). 

(e) An approximate number of particles that provide the desired porosity 
is calculated (section 2.5). 

(f) The radii of the particles are samples of independent identically dis- 
tributed random variables obeying the adopted distribution (section 2.6). 

2. The particles are placed by a Simple Sequential Inhibition - SSI process on 
the arbitrary closed region (section 2.7). The number of particles is deter- 
mined iteratively as to produce an approximation of the desired porosity. 

The result of these steps is a plausible model of particle distribution in 
granular media, obtained with low computational cost due to its geometric 
nature. Algorithm 1 presents, in pseudocode, an efficient implementation of 
this procedure. We show that this procedure is able to simulate soils of low 
porosity, which are usually obtained by dynamical techniques that are much 
more time consuming. Such models can be used as input data for numerical 
simulations of problems with discontinuous media using the discrete element 
method. 

The SSI process we employ can be seen as a RSA process with the following 
characteristics: the particles only interact by exclusion in a compact planar 
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domain, i.e., there is neither attraction nor repulsion, no order is introduced 
in the process, i.e., the final state is jammed, particles arrive sequentially and 
there is no relaxation, i.e., the (fc + l)th particle enters the sample only if the 
fcth particle does. 

This paper unfolds as follows. Section 2 presents the methodology of each 
stage of our proposal, including pseudocode algorithms. Section 3 discusses 
the main results, while section 4 presents the conclusions and future venues of 
research. 



2 Particle size analysis and generation 

2.1 Granulometric curve and porosity 

A granulometric curve measures the percentage of the sample that falls into 
pre-established ranges of grain sizes. This information is obtained by seiving 
and hydrometer analyses, for instance. The results of such measures is provided 
in tabular form: (d^, Ci)i<i<£), being q the (cumulative) proportion of particles 
whose diameter is less than di, and D the number of diameters considered. 

The porosity r] G [0, 1] is the ratio between the empty space in the sample 
to the total sample volume. It is estimated through the determination of the 
bulk density of the porous sample, then determining the density of the skeletal 
material, and finally correlating density and volume. 

The granulometric curve and the porosity are the only required input for 
the rest of the methodology. 



2.2 Histograms 

The granulometric curve is used to form a tractable and expressive histogram 
of particle sizes. As presented in Figure 2, the diameters are not evenly spaced 
and, thus, the intervals are not equal. Such unevenness is alleviated in the 
semilogarithmic scale, so our proposal uses such transformation. 

Consider the original data (rfi,Ci)i<i<£) and its representation in semilog- 
arithmic scale {£i,Ci)i<i<D, where ii = \ogdi for every 1 < i < D. The log- 
histogram of the data can be formed with the midpoints of pairs of contiguous 
log-diameters and the corresponding proportion of particles, i.e., {mi,pi)i<i<D-i, 
where nii — {ii+i +£i)/2 and pi = q+i — a. These log- histograms will be used 
to obtain more data from pseudorandom sampling, in order to obtain para- 
metric models for the data. 

Each dataset is, therefore, transformed into a log-histogram. Instead of us- 
ing the pairs (mi,pi)i<i<D_i, we draw Ni independent identically distributed 
samples from the Uniform distribution on with = [kpi] and k a 

convenient number, for instance k ^ 1000. With this, we end with X^tTi^ — 
M pseudorandom diameters d — (di, . . . ,dM) which, if given as input, would 
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produce the original granulometric curve {di,Ci)i<i<D- In this way, an arbi- 
trary number of samples can be obtained from a single granulometric curve, 
leading to very precise parameter estimation. 

These data d will be then used for computing estimates of parametric 
models able to describe the original data. 



2.3 Models 

2.3.1 Hyperbolic distribution 

The unidimensional hyperbolic distribution is a continuous probability law 
characterized by the fact that the logarithm of the probability density function 
is a hyperbola. It was introduced by [5] as a model for the log-size distribution 
of sand, based on the studies of [3] about aeolian sand deposits. 

This distribution and its multidimensional extensions have many appli- 
cations in a variety of fields: geology, astronomy, fluid mechanics and eco- 
nomics [4,6,7,9], to name a few. 

One of the representations of the probability density function, presented 
in [8], is: 

v/a2 - /32 exp{~a^S^ + {x - + /3{x - ^)} 
f(x;a, p, u,d) = , , 

where x,^,l3 e R, S,a G R+ such that < |/3| < a, and Ki, i = 1,2,3, 
denotes modified Bessel functions of the third kind. The location parameter 
fj, is the abscissa of the point of intersection between the linear asymptotes 
of the hyperbola, while 6 is the scale parameter. The two parameters a and 
/3 determine the shape, being a resposible for the steepness and l3 for the 
skewness. The distribution is symmetric for f3 = 0. 

If is a hyperbolic distributed variable, its expected value and mode are 

EiX) ^ , + , SpK^jSy^) ^ ^ 



and its variance is 



Var(X) ^SK^iSV^^^^ 



Another representation of the density function of the hyperbolic distribution 
is: 

exp{-c(Vr + ^V'r7(^-7r^)} 
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where a;, /i, tt € R, and (5, C G R+. The parameter C is a measure of the degree 
of peakedness, and tt is a measure of the asymmetry of the distribution. When 
TT — 0, the distribution is symmetric. Together, they determine the shape 
of the distribution. ^ and d are parameters of location and scale. This last 
parametrization is employed by the HyperbolicDist package of R [39], which 
also provides techniques for computing maximum likelihood estimates. 



2.3.2 Gamma distribution 

The gamma distribution is a continuous probability distribution of two pa- 
rameters. It has been used in a wide range of disciplines: from climatology (in 
particular, to study rainfall; see [1]), to material analysis [10], to name a few. 

It has a scale parameter A and a shape parameter a, both positive, in its 
density function: 

where F is the incomplete gamma function given by r{a) ~ x°'^^e^^dx. 
Alternatively, it can be parameterized in terms of the rate f3 — 1/A. 

If X is a gamma distributed variable, its expected value, mode and variance 

are: 

E{X) = aX, Modc(X) = (a - 1)A, and Var(JSi:) = aA^ 
respectively. 



2. 3. 3 Lognormal distribution 

The lognormal distribution is the distribution of any random variable whose 
logarithm is normally distributed. It is widely used in physics, statistics, geol- 
ogy, economics, biology etc. Its probability density function is 

where /i G R and a > 0. 

If X is a lognormally distributed variable, its expected value, mode and 
variance are 

E{X) = e^+"'/2^ Mode(X) = e^-"', and Var(X) = (e"' - l)e^''+^\ 
respectively. 
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2.3.4 Weibull distribution 

This law is often called the Rosin-Rammler distribution when used to describe 
the size distribution of particles [40]. The probability density function of a 
Weibull random variable X is 

f{x- a, A) = aA"x"-ie-(^")° (x), 

where a, A > are the shape and the scale parameters, respectively. 

If X is a Weibull distributed variable, its expected value, mode and variance 
are: ^ 

E{X)^Xr(l + ^, Mode(X) ^ a(^— ^)°, 

and 

respectively. 



2.4 Choice of the best model 

Each sample is used to estimate the maximum likelihood parameters 9 that 
index the four distributions 2? considered as candidate models, and then the 

goodness-of-fit test is applied. In principle, the model with highest p-value 
should be chosen as the best statistical description and, in case of tiers, the 
one with least computational requirements should be the final choice. 

Other goodness-of-fit test as, for instance, the Kolmogorov-Smirnov test, 
can be used, as well as other estimation procedures, e.g., moment or analogy 
estimators. 



2.5 Number of particles 

We will derive an initial guess for the number of particles to be generated in 
order to have a good approximation of the desired prosify. 

Assume the closed area of interest, W, has size n{W), which we want to 
fill with particles whose radii follow the parametric model 2?(6') chosen as 
previously discussed. We also have the desired porosity 77, estimated from the 
soil sample: 

^ ^ li{W) - Vp 

where Ve and Vp are the empty and particles areas, respectively. The area of a 
disk with radius R, distributed according to the law 15(0) is a random variable 
given hy V = ttR^. Then, assuming that N particles are placed within W, 

N 

Vp = 7rJ2Rl 

i=l 
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Using expected values, 

N N 
i=l 1=1 

and both variance and expected value are given explicitely as functions of the 
estimated parameter 0, so 

^ ^ {i-vU W) 



This number of particles N is intended as a rough estimate of the needed 
number of radii required to fill W with the desired porosity and granulometry. 



2.6 Generation of particles 

Four distributions are considered in this work, namely the Hyperbolic, Gamma, 
Log-normal and WeibuU laws. Sampling from the two first requires specialized 
algorithms, samples from the third are the exponential transformation of nor- 
mal deviates, and inversion is required in order to sample from the last one. All 
these samples are produced calling functions available in the R platform [17]. 

Algorithm 1 presents, in pseudocode, a sequential approach for obtaining 
a good approximation to the desired porosity 77, while keeping the statistical 
properties of the particles. In this algorithm, if > 1 is a factor which controls 
the number of extra particles to be generated; we set K — 10, a. good compro- 
mise between the relatively low cost of generating outcomes versus the high 
cost of initializing the generator. 

This algorithm places a number of particles obeying the distribution 
which approximates the desired porosity i] on the region of interest W. 

It is noteworthy that the presented methodology can be immediately ex- 
tended to compact volumes. 



2.7 Spatial distribution by the SSI process 

Line 6 of Algorithm 1 makes a call to the procedure presented in Algorithm 2, 
which implements the sampling from the simple sequential inhibition point 
processes. Such processes are defined in terms of a window W, the number of 
points and the exclusion radii among them [11]; in our case, the number of 
points is, instead of specified a priori, controlled by the desired porosity 77 on 
W. 

The Simple Sequential Inhibition (SSI) point process for n points on W 
with exclusion radii r = (ri, . . . , r„) and maximum number of iterations jmax 
tries to place n non-overlapping disks of radii r by testing sequentially each 
disk against the previous ones. The first disk is placed uniformly on W, pro- 
vided it does not surpass the boudaries. Subsequent steps sample a point in W 



Stochastic particle packing with specified granulometry and porosity 



11 



Algorithm 1 Sequential sampling of radii for a desired porosity 



1: procedure SEQUENTlALRADll(Ty , A'", r], V(9), K) > Region of interest, number of 

particles, porosity, model and factor 
2: Obtain r = (ri, . . . ,rKN) samples from independent identically distributed random 

variables following the "Did) distribution 
3: Allocate the vector of coordinates c = {{x\, yi), . . . , {xxn ^ Vkn)) 
4: i 0, »7i <— 1 > Initialize counter and porosity 

5: while rn > ri do 

6: condition <— SSI(iy, i, r, c) > Defined on Algorithm 2 

7: if condition = TRUE then 

8: i <— i + 1 > Update counter 

9: Tji Tji-i — Trr^ / ^{W) > Update porosity 

10: else 

11: break o No more particles will be placed 

12: end if 

13: end while 
14: if j = then 

15: return FAIL > The first particle could not be placed in W 

16: else 

17: return yi , . . • , {xi,yi,ri),r)i) l> Returns coordinates, radii and current 

porosity 
18: end if 
19: end procedure 



uniformly and independently of previous disks, and verifies if there is no over- 
lapping; if there is not, the disk is placed, otherwise up to jmax independent 
trials are made. The algorithm stops when all the disks have been placed, or 
when the maximum number of iterations has been reached, whichever takes 
place first. Our implementation of the SSI point process is presented in Algo- 
rithm 2. 



Algorithm 2 Simple sequential inhibition point process 

procedure SS\{W,i,r , c) t> Region of interest, iteration, radii and coordinates 

Define j'max > The maximum number of trials per particle 

while j < Jmax 

do 

Sample {x, y) uniformly on W 
if mini<fc<i{d((x,iy),(xfe,yfe)) > rj then 

c(i) <— {x,y) > Place particle i at {x,y) 

return TRUE l> Success 

else 

i ^ J + 1 

end if 
end while 

return FALSE > Unable to place the i-th particle 

end procedure 



Notice that the order in which the particles are places is not altered. In 
particular, if particle i < 2 does not fit after jmax trials, the algorithm returns 
the current state which consists of the previous i — 1 particles. Even if particle 
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i+1 fits trivially the current state, it is not included. Selecting particles by their 
size would imply independent identically distributed random variables model, 
and switching to a collection of correlated deviates: the order statistics. 

Dense samples are, typically, harder to obtain than situations of high poros- 
ity. In those cases it is likely that the denser the desired outcome the more 
results will have to be discarded until a sample with the desired property is 
obtained. 

3 Results 

This section presents results related to the two main steps already discussed, 
namely, particle size generation and placement. The example used to illustrate 
this work is based on collected disturbed samples from the Barreiras formation 
and other fluvial deposits in Alagoas state, Brazil, in order to study their 
granulomctric properties. 

The granulomctric information is obtained in tabular form: {di, Ci)i<i<D, 
being q the (cumulative) proportion of particles whose diameter is less than 
di, and D the number of diameters considered. 

The samples used as input for this analysis are presented in Table 1. Fig- 
ure 2 presents the result of the sieving process applied to Sample 1. The table 
shows the diameters, in decreasing order and in mm, and the corresponding 
cumulative percentage of passing particles. This same information is shown 
in the form of a granulomctric curve in linear scale (right top); the sieves 
diameters are show superimposed to the abscissa axis. The clutter effect of 
small sieves is evident, therefore the frequent choice of presenting such curves 
in semilogarithmic scale (right bottom). Notice that the diameters are not 
evenly distributed deserving, thus, a careful treatment. 

Table 1 Soil types analized 



Sample 


Soil type 


Porosity rj 


1 


sand, silt with traces of de gravel 


0.35 


2 


sand with traces of silt 


0.43 


3 


sand, portions of clay, traces of silt and gravel 


0.42 


4 


sand with traces of silt 


0.39 


5 


sandy clay with portions of silt 


0.40 


6 


sandy clay with portions of silt and traces of gravel 


0.33 


7 


sand 


0.44 



Each data set from Table 1 was analyzed using the procedure presented 
in previous sections, i.e., (i) its granulomctric curve was obtained by sieving 
(see Figure 2 as an example), (ii) a histogram was formed (Figures 3 and 4), 
(iii) samples of pseudorandom diameters were drawn. These diameters were 
then analyzed by fitting the gamma, lognormal, Weibull and hyperbolic distri- 
butions. Maximum likelihood estimation was performed using the f itdistr 
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function from the MASS package of R (for the gamma, lognormal and Weibull 
distribution) and the hyperbFit function from the HyperbolicDist package 
(for the hyperbohc distribution). The estimated parameters, rounded to the 
third decimal place, are given in Table 2. 



Table 2 Estimated parameters 



Gamma Lognormal Weibull Hyperbolic 



cii 

CO 


a 


/3 




fj 


a 


A 


77 


c 


S 


M 


1 


43.717 


4.170 


2.338 


0.159 


10.281 


11.046 


-1.055 


0.750 


0.403 


11.892 


2 


342.999 


32.155 


2.366 


0.054 


18.806 


10.943 


0.235 


15.073 


2.072 


10.132 


3 


81.896 


7.437 


2.393 


0.113 


11.241 


11.513 


-3.369 


19.473 


1.367 


15.977 


4 


387.377 


36.738 


2.354 


0.051 


21.186 


10.799 


-0.138 


705.286 


14.081 


12.495 


5 


114.429 


11.722 


2.274 


0.095 


12.564 


10.163 


-1.171 


45.259 


3.816 


14.381 


6 


90.818 


9.153 


2.289 


0.106 


10.231 


10.381 


0.107 


2.635 


1.304 


9.698 


7 


334.601 


130.210 


2.403 


0.055 


16.146 


11.372 


0.159 


2.025 


0.631 


10.895 



The p-values of the test were computed, and they all resulted above 
0.99 with one exception, the Weibull distribution for sample 7. Therefore, for 
the samples here analyzed, the choice of the "best" model can be guided by 
the computational cost of producing pseudorandom deviates or by any other 
criterion. 

We simulated samples from each of the types of soil previously analized 
using the models presented in Table 3, which describe the logarithm of the 
diameters. These simulations were performed in squared boxes specifying the 
desired porosity. Table 3 also presents the number of particles generated in 
each case {N), along with the desired and obtained porosities (77 and rj, re- 
spectively). Given the similarity of the results, in Figure 5 we only present the 
simulation corresponding to Sample 6. Figure 5(a) presents the 34688 particles 
that were generated to fill a box with the specified porosity, while Figure 5(b) 
shows a detail where both the varibility of the model and the non-overlapping 
results are enhanced. 



Table 3 Chosen model for each sample, number of particles, desired and obtained porosities 



Sample 


V 


N 


V 


V 


1 


WeibuU 


6704 


0.350 


0.345 


2 


Hyperbolic 


4846 


0.430 


0.421 


3 


Gamma 


1922 


0.420 


0.418 


4 


Lognormal 


7840 


0.390 


0.390 


5 


WeibuU 


27486 


0.400 


0.399 


6 


WeibuU 


34688 


0.330 


0.328 


7 


Hyperbolic 


1611 


0.440 


0.430 
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4 Conclusions and future work 

We proposed and assessed a fast geometric algorithm which produces spatially 
homogeneous samples with specified porosity and granulometric curve using 
disks. Starting with the granulometric curve and porosity of a soil, our algo- 
rithm produces samples of any size in closed arbitrary domains. The outcomes 
of this algorithm can be used as granular media samples in computational 
simulations, e.g., discrete elements. 

The algorithm allows the use of any suitable distribution for the particle 
size. It was assessed with seven samples and four distributions, and in every 
case it produced samples wich mimic well the input data. 

The proposed strategy can also be used in the simulation of other physical 
systems. Kadau et al. [30] propose a twodimensional contact dynamics model 
as a microscopic description of a collapsing suspension/soil to capture the es- 
sential physical processes underlying the dynamics of generation and collapse 
of the system. Kadau and Herrmann [29] study the influence of the granular 
Bond number on the density profiles and the generation process of packings, 
generated by ballistic deposition under gravity. Both works deal with loose 
granular media and, thus, there is no need to impose strong physical con- 
straints. 

The research continues with the proposal of a hybrid technique which, 
starting with our geometrical algorithm, produces a final configuration in geo- 
static equilibrium by relaxation in both 2D and 3D domains. 
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Diameter [mm] 


Cumulative % 


4.800 
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98.65 


1.200 


83.81 


0.840 


70.66 


0.710 


61.49 
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56.29 
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0.420 
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36.41 
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20.49 
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0.025 
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1.09 


0.002 


0.00 



Cumulative proportion 



20 40 60 80 100 




(a) Original granulometric curve 
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(b) Granulometric curve in semilogarithmic scale 

Fig. 2 Granulometric data and curves: linear (a, top) and semilogarithmic (b, bottom), 
sample 1 



18 



A.C. Prery et al. 




Fig. 3 Loghistogram and fitted density for sample 1 
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(e) Sample 6 (f) Sample 7 

Fig. 4 Loghistograms and fitted densities 
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